#################################################################################
#################################################################################

###          R code to replicate simulation experiment reported in            ### 
###      Supplementary materials (Appendix E; Figure E1: Panels A-C) to       ###
###  "The Index of Emancipative Values: Measurement Model Misspecifications"  ###

#################################################################################
#################################################################################


### Load required R packages ###

library(simsem)
library(lavaan)
library(semPlot)


### Generate true model ### 

genmod <- "
f1 =~ 1*y1 + 1*y2 + 1*y3  
f2 =~ 1*y4 + 1*y5 + 1*y6 
f3 =~ 1*y7 + 1*y8 + 1*y9  
outcome  ~ .7*f1 + .3*f2 + (-.3)*f3  
outcome  ~~ 1*outcome 
f1 ~~ 0.6*f2
f1 ~~ 0.6*f3
f2 ~~ 0.6*f3
"


### Simulate data ###
set.seed(666)
dat <- simulateData(genmod, sample.nobs = 100000)



### Estimate true CFA model ### 

# define model
mod <- "
f1 =~ 1*y1 + NA*y2 + NA*y3 
f2 =~ 1*y4 + NA*y5 + NA*y6 
f3 =~ 1*y7 + NA*y8 + NA*y9  
outcome ~ NA*f1 + NA*f2 + NA*f3 
"

# estimate model
out <- cfa(mod, dat)

# return results 
summary(out, standardized = TRUE, fit.measures=TRUE, rsq=TRUE)



### Model 1 Misspecified (Reflective, with Second-order factor) ### 

mod2 <- "
f1 =~ 1*y1 + NA*y2 + NA*y3 
f2 =~ 1*y4 + NA*y5 + NA*y6 
f3 =~ 1*y7 + NA*y8 + NA*y9  
f4 =~ 1*f1 + NA*f2 + NA*f3  
outcome ~ NA*f4 
"

out2 <- cfa(mod2, dat) 
summary(out2, standardized = TRUE, fit.measures=TRUE, rsq=TRUE)


### Model 2 Misspecified (Formative Index) ### 

# Create index
dat$y.agg <- (dat$y1 + dat$y2 + dat$y3 + dat$y4 +dat$y5 + 
dat$y6 + dat$y7 + dat$y8 + dat$y9)/9

mod3 <- "
outcome ~ y.agg 
"

out3 <- cfa(mod3, dat) 
summary(out3, standardized = TRUE, fit.measures=TRUE, rsq=TRUE)



### Figure E1: Panel A ### 

tiff('Figure E1a.tiff', units = "cm", width = 16, height = 12, res = 600)

semPaths(out, "est", 
layout = "tree", residuals = FALSE,
edge.label.cex= 1, 
thresholds = FALSE, 
intercepts = FALSE, 
nodeLabels= c("Y1", "Y2", "Y3", "Y4", "Y5", 
"Y6", "Y7", "Y8", "Y9", "Outcome", "F 1", "F 2", "F 3"),
optimizeLatRes = TRUE,
sizeMan = 9,
sizeMan2 = 5,
sizeLat = 7,
mar = c(1,2,1,2),
label.cex= 0.75, 
label.scale= FALSE, 
esize= 2,  
weighted= FALSE,
edge.color= "black") 

dev.off()


### Figure E1: Panel B ### 

tiff('Figure E1b.tiff', units = "cm", width = 16, height = 13, res = 600)

semPaths(out2, "est",
layout = "spring", residuals = FALSE,
edge.label.cex= 1, 
thresholds = FALSE, 
intercepts = FALSE, 
nodeLabels= c("Y1", "Y2", "Y3", "Y4", "Y5", 
"Y6", "Y7", "Y8", "Y9", "Outcome", "F 1", "F 2", "F 3", "F 4"),
optimizeLatRes = TRUE,
sizeMan = 9,
sizeMan2 = 5,
sizeLat = 7,
mar = c(1,2,1,2),
label.cex= 0.75, 
label.scale= FALSE, 
esize= 2,  
weighted= FALSE,
edge.color= "black")

dev.off()



### Figure E1: Panel C ### 

tiff('Figure E1c.tiff', units = "cm", width = 10, height = 8, res = 600)

semPaths(out3, "est", residuals = FALSE,
edge.label.cex= 1.75, 
thresholds = FALSE, 
intercepts = FALSE, 
nodeLabels= c("Outcome", "Composite\nScore"),
optimizeLatRes = TRUE,
sizeMan = 20,
sizeMan2 = 12,
label.cex= .8, 
mar = c(5,5,5,5),
label.scale= FALSE, 
weighted = FALSE,
esize= 3.5,  
border.width = 2.5,
edge.color= "black")  

dev.off()

# rm(list = ls())
# q()
